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Langley Stability and Transition Analysis Code (LASTRAC) is a general-purpose, physics-based transition 
prediction code released by NASA for laminar flow control studies and transition research. This paper 
describes the LASTRAC extension to general three-dimensional (3D) boundary layers such as finite swept 
wings, cones, or bodies at an angle of attack. The stability problem is formulated by using a body-fitted 
nonorthogonal curvilinear coordinate system constructed on the body surface. The nonorthogonal coordinate 
system offers a variety of marching paths and spanwise waveforms. In the extreme case of an infinite swept 
wing boundary layer, marching with a nonorthogonal coordinate produces identical solutions to those 
obtained with an orthogonal coordinate system using the earlier release of LASTRAC. Several methods to 
formulate the 3D parabolized stability equations (PSE) are discussed. A surface-marching procedure akin to 
that for 3D boundary layer equations may be used to solve the 3D parabolized disturbance equations. On the 
other hand, the local line-marching PSE method, formulated as an easy extension from its 2D counterpart 
and capable of handling the spanwise mean flow and disturbance variation, offers an alternative. A linear 
stability theory or parabolized stability equations based N-factor analysis carried out along the streamline 
direction with a fixed wavelength and downstream-varying spanwise direction constitutes an efficient 
engineering approach to study instability wave evolution in a 3D boundary layer. The surface-marching PSE 
method enables a consistent treatment of the disturbance evolution along both streamwise and spanwise 
directions but requires more stringent initial conditions. Both PSE methods and the traditional LST 
approach are implemented in the LASTRAC.3d code. Several test cases for tapered or finite swept wings and 
cones at an angle of attack are discussed. 


I. Introduction 

I nstability of laminar flow and transition to turbulence in a three-dimensional boundary layer such as a finite-span 
swept wing or axisymmetric forebodies at an angle of attack are important issues of modem aerodynamic designs. 
At cruise conditions, the drag caused by turbulent boundary layers over the surface accounts for a large portion of 
the total drag. Proper laminar flow control (LFC) to prevent or delay the onset of transition would significantly 
increase the vehicle performance and result in cost savings. The recent Defense Advanced Research Projects 
Agency (DARPA) effort on the Quiet Supersonic Platform (QSP) is aimed at developing technologies to achieve a 
low-boom, laminar flow design at supersonic speed [1]. In conjunction with the NASA and industry Supersonic 
Vehicle Technology (SVT) project, these transition-related efforts are geared towards the development of enabling 
technologies that will eventually lead to the next generation supersonic vehicles beyond Concorde. The core effort 
of these projects involves computational tool development and wind tunnel or flight experiments. 

At NASA Langley Research Center, recent transition-related works have been focused on the development of a 
physics-based, integrated transition prediction framework that involves theoretical and numerical methods, 
prediction software development, and experimental investigations. This integrated tool consists of software that is 
designed to handle various stages of the transition process. The software encompasses a range of methods from the 
conventional N-factor correlation (based on either the linear stability theory (LST) or the more advanced parabolized 
stability equations (PSE) approach) to simulations based on the absolute disturbance amplitudes. The Langley 
Stability and Transition Analysis (LASTRAC) software, made publicly available in 2002[2, 23]. represents the first 
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release of the integrated tool development. The code can be used for transition-related efforts for 2D or infinite 
swept wing configurations. This paper is focused on a major extension of LASTRAC to general 3D boundary layers. 

Stability of 3D boundary layers has drawn much attention for several decades. Most of these investigations were 
based on the linear stability theory (e.g. Nayfeh [3], Cebeci & Stewartson [4], and Malik & Balakumar [5]) where 
stability analysis is done locally by assuming a slow-varying mean flow on the 3D surface. An eigenvalue problem 
similar to that in 2D boundary layers was formed in these analyses. N-factor integration was then carried out to 
represent the disturbance growth. The computed N-factor distribution on the 3D surface is thus closely related to the 
integration path chosen. Either the inviscid streamline or the group velocity line was commonly chosen to integrate 
the disturbance growth. Most of transition correlation results did not favor one over the other (see Malik & 
Balakumar [5]). More recently. Atobe and Chen [6] also used a similar method for transition prediction on a Mach 2 
natural laminar flow wing model. 

The PSE method has been proven to be a viable alternative to traditional linear theory for 2D and infinite swept 
wing boundary layers. The PSE approach offers a very efficient marching solution that can account for the effects of 
boundary layer growth and nonlinear mode interaction. The capability of handling the latter effect forms the basis of 
a much more efficient method for absolute amplitude transition prediction without resorting to direct numerical 
simulations (DNS). The extension of PSE to general 3D boundary layers, however, is nontrivial due to several 
reasons. Firstly, unlike in the infinite swept wing boundary layers, the mean flow varies along the spanwise direction 
and this variation prevents the disturbance from being expressed as a simple Fourier decomposition in the spanwise 
direction. Therefore, a finite-difference description in physical space or a Fourier collocation scheme must be used 
to incorporate the changes along the spanwise direction. Secondly, proper disturbance boundary conditions are not 
easily identified in the spanwise direction except in the cases where periodic or symmetry conditions may be used 
(such as in an axisymmetric cone). Thirdly, because of the added degree of freedom on a 3D surface, the 
“stream wise 1 ' and “spanwise” directions are less apparent. As a result, a proper waveform is more difficult to 
construct as than its counterpart in a 2D boundary layer. 

To circumvent the above difficulties, several strategies with varying degree of complexity may be employed. With 
proper disturbance boundary conditions, a plane-marching approach may be the most intuitive to use along a pre- 
identified streamwise direction (e.g., along the free stream). In contrast to the 2D cases, the modal shape function is 
now a two-dimensional function. Given an initial distribution on a cross plane, the marching solution may be carried 
out. In the case where disturbance boundary conditions are difficult to find, a window or fringe function approach as 
suggested by Spalart and Watmuff [7] and Guo et al. [8] may be used to render a finite spanwise domain periodic 
such that a Fourier decomposition may be used. It should be noted that even with proper spanwise boundary 
treatment, additional issues such as the determination of a proper streamwise wave number also arise. Compared 
with 2D cases, this new approach is not as efficient because a cross-plane solution would require solving a large- 
bandwidth matrix. Domain decomposition in conjunction with parallel processing or disk caching techniques are 
necessary to achieve more acceptable numerical efficiency. 

For efficiency reasons, it is more desirable to have a simplified “local” approach that does not require a plane 
solution. Herbert [9] suggested a streamwise and spanwise surface-marching approach. In this method, two complex 
wave numbers are defined in accordance with two pre-defined spatial coordinates on a 3D surface. The evolution of 
the shape function is solved iteratively with two normalization conditions to determine two wave numbers. Initial 
conditions along two initial lines are required to start the marching process. The distribution of the disturbance 
growth may be obtained by integrating both wave numbers simultaneously on the surface and accounting for 
additional contributions coming from the shape function evolution. Because both wave numbers vary during the 
marching process, there is no clear definition of a particular disturbance mode. It appears that a growth pattern may 
be strongly related to the initial condition provided. 

Alternatively, a local PSE approach may be constructed by locally varying the spanwise coordinate to align with a 
direction where the mean flow variation is minimal. Following the lead from the LST approach, a PSE marching 
solution along the inviscid streamline with locally varying spanwise direction offers a more economical method to 
account for 3D effects. In this paper, we formulate the stability problem by using a body-fitted nonorthogonal 
coordinate system that allows a flexible variation of spanwise coordinates along the marching direction. As will be 
discussed in the following section, mean flow and disturbance variation along the variable spanwise direction may 
be accounted for in this approach. In many applications, this spanwise variation does not play a significant role on 
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the overall growth and may be neglected. In cases where substantial spanwise mean flow variation is present, the 
effect can be accounted for in this local analysis. 

Both the surface-marching and local line-marching PSE methods for 3D boundary layers are investigated in the 
present study. These 3D PSE methods along with the more traditional LST-based approach are implemented into the 
new release of the LASTRAC.3d code. These different methods provide a means for cross validation. One of the 
main objectives of this paper is to justify the use of a local line-marching PSE approach for instability wave 
evolution in a 3D boundary layer. Because of its numerical efficiency, this local approach is recommended as an 
engineering tool for transition correlation. The flexible nonorthogonal coordinate system allows a varying spanwise 
direction that leads to better transition correlations. In this paper, we discuss the problem formulation for several 3D 
PSE methods. Some preliminary linear results for a transonic tapered wing, a finite 3D swept wing, and an 
axisymmetric cone at an angle of attack are included. Only linear results are presented in this paper. Transition 
prediction issues are assessed based on the disturbance amplification (N-factors). 


II. Problem Formulation 

To allow a more general waveform for 3D boundary layers, we use a nonorthogonal curvilinear coordinate system 
(x 1 , X 2 , X 3 ) where X and x 3 are along the boundary layer surface and X 2 along the wall-normal direction (see 

Fig. 1). Coordinates X 1 and x 3 are roughly aligned with the streamwise and spanwise directions, respectively and 
they are nonorthogonal in general. The Cartesian coordinate system is denoted by (x, V, -) . In the (x',x 2 ,x 3 ) 
coordinate, the distance between two points infinitesimally separated in space may be evaluated by 

ds 2 = g iJ dx'dx J 

where g t] is the fundamental metric tensor. The determinant of the matrix with an element of g tj is equal to the 
square of the Jacobian of the coordinate transformation. The contravariant metric tensor is defined according to 

g»g* = 8 i 

where 8- is the Kronecker delta and the tensor summation is implied. Both g tj and g' r can be evaluated based on 
the coordinate transformation from (x,y,z) to (x*,X ,X 3 ). 


The Navier-Stokes equations in the curvilinear system can be written as (Aris [10]) 

dp 


dt 


Hpv’h= 0 


do* 

dt 
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where p, p,T , p,A,K and C p are density, pressure, temperature, first and second coefficient of viscosity, 

thermal conductivity and constant pressure specific heat, respectively. The quantity v' is the contravariant velocity. 
The divergence of any vector is defined as 
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where the Christoffel symbol of the second kind is indicated by the curly bracket. The second-order tensor l/ ; is 
defined by 


V 


m j ' L* 

dxj l^yj 


The viscous dissipation function may be evaluated by 



+ g""o 


J \g 

>m f\c> ir 


V 



For a Cartesian coordinate system, the metric tensor g t] reduces to 

g =5' 

* y j 

For an orthogonal curvilinear coordinate system, the tensor becomes 

= h 2 Sj (no summation over /). 

The quantity h t is the metric associated with the curvature along the x‘ direction. In either case, g lf becomes a 

diagonal tensor and the governing equations reduce to the conventional form. In general, g t/ is a full (but 
symmetric) second-order tensor. 

We decompose the dependent variables <p t (p,v\v 2 ,U 3 ,7’) into a steady-state mean part plus a perturbation, 
(p t {x x ,x 2 ,x 3 ,t) - ,x 2 ,x i ) + </>{x x ,x 2 ,x 3 ,t) (2) 

If we substitute the above equation into Eq. (1) and subtract the steady-state equations, then we have the governing 
equations for the disturbances as 


dt 



+ B^- + C^- + D</> = V XX 
dy oz 


d 2 <t> 
dx 2 + 


v ?_ U v **- 

” dx dy " dy 1 y ~ dydz dxdz 



( 3 ) 


Here, for convenience, we have replaced the coordinates (x',A' 2 ,X 3 ) with (x,y,z) . The coefficient matrices 
r, A,B,C , etc. are functions of all three coordinates for 3D boundary layers. They consist of a steady-state 
(functions of mean quantities) and time-dependent parts (functions of unsteady perturbations), e.g.. 
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r=r 'w+r^, $) 

The perturbation can be expressed in the following discrete Fourier series with a fundamental disturbance frequency 
(O 


M 


</>(x,y,z,t) = £ y/ m {x,y,z)e~ 

m~-M 

The governing equation for a single Fourier mode becomes 


( 4 ) 


A^B^C^HD-imaryr.^V^ + V + 

dxdydz dx 2 x> dxdy 

V n ^ + V^ + V^ + V_^ + F , ( 5 ) 

■ ay oyoz oxoz dz 

where the forcing function F m contains contributions from other Fourier modes due to nonlinear interaction and can 
be neglected for small disturbances. The coefficient matrices contain only the steady-state part ( r ? , A 1 , etc. ) and the 
time-dependent part ( T", A” . etc.) has been moved to the right-hand side forcing function. 


For a 2D or an infinite swept wing boundary layer, the basic flow is invariant along the spanwise direction. Periodic 
boundary condition in z is thus assumed. The mode shape function may be further decomposed into a discrete 
Fourier series. 


¥, 


,(*,T,5)= £ V' nm (x,y)e m/,: 


n=-N 


( 6 ) 


Substitute Eq. (6) into Eq. (5), we have a set of coupled partial differential equations (PDE) for y/ mn in two space 

dimensions, (x, y). The derived equation is the harmonic linear Navier-Stokes (LNS) equation solved by Streett [11] 
and Guo et al. [12] using a large-matrix solver. Further assumptions based on scaling or parabolizing the derived 
equations may be made to simplify the solution procedure. In the PSE approach, a streamwise wave number is 
introduced to minimize the parabolizing error and to reduce the stringent resolution requirement in x. These 
procedures are discussed in the previous release of LASTRAC in details (Chang [2, 23]). 

Two major issues stand in the way if we attempt to apply Eq. (6) directly for a full 3D boundary layer. Firstly, a 
periodic condition in r is clearly invalid. In fact, proper boundary conditions are often difficult to identify in the 
spanwise direction. Depending upon the coordinate system used, there exists no clear physical distinction between x 
and z in some cases. For the ease of the following discussion, we assume that the coordinate x is roughly aligned 
with the free stream direction and the spanwise coordinate (which may or may not be perpendicular to x in general) 
may be chosen to align with a direction where the mean flow variation is minimal. The second issue is that for a 
given spanwise direction, the wave number p is in general complex for a full 3D boundary layer. The mean flow 
variation in z allows a disturbance mode to grow or decay in that direction. To overcome the above two issues, 
several methods with a varying degree of complexity may be used to solve Eq. (5) in a 3D space. 

A. Plane-Marching PSE 

In some cases, proper boundary conditions do exist in the spanwise direction. For example, periodic or symmetry 
boundary conditions may be applied in the azimuthal direction for axisymmetric configurations such as cones. 
Dirichlet or Neumann boundary conditions and no-slip conditions may be applied to the free stream and wall, 
respectively. 
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Similar to the 2D cases, with a quasi-parallel assumption along the streamwise ( x ) directions, we can rewrite y/ m as 


¥ m = v m {y^)^ 


( 7 ) 


where the streamwise wavenumber CC m is a function of the disturbance frequency according to the dispersion 
relation a m = OL m (mco) . The linearized governing equation becomes 


dy 



+ Dy/ m 


= y — ^- + F 

r v\ -V 2 yz 

- dy 


d 2 w d 2 w 

r m _j_ y ^ Y m 

dydz " 8z 2 


( 8 ) 


With homogeneous boundary conditions on the y-z plane discussed above, a 2D eigenvalue problem similar to that 
investigated for secondary instability of stationary crossflow instability by Malik et al. [13] may be formed with 

y/ m (y,z) as the 2D eigenfunction. Solving a 2D eigenvalue problem is rather time-consuming even with a 

moderate grid size on they-; plane. Numerical efficiency may be significantly improved by using Amoldi’s method 
as was done by Lin and Malik [14] for the instability of the attachment line boundary layer. 


To further account for the nonparallel effect, we may introduce a streamwise varying wavenumber, CC m (x ) , such 
that 


~ if a(£)d£ 

V' m =V' m (x,y,z)e « 


( 9 ) 


The governing equations become 

dx dy 8z " dy>~ ' dydz 8z 


( 10 ) 


where additional terms associated with the shape function derivatives in x are introduced as compared to Eq. (8). 
Similar to its counterpart in 2D cases, Eq. (10) may be solved by a plane-marching procedure in X provided that 

initial conditions are given on a cross plane. Unlike in 2D cases, determining the streamwise wave number a m 
remains an issue. The introduction of the streamwise wave number a m in the PSE approach is, primarily, to resolve 
the oscillatory (fast scale) part of the instability wave. Therefore, a constant a m for the whole (y, z) plane becomes 
inadequate when the disturbance mode structure varies substantially on the plane. We can make a m a function of x 

and z. The distribution of the streamwise wave number must then be determined by a normalization procedure local 
to the spanwise location r. 

For flow configurations such as a finite 3D wing where periodic conditions are not applicable, a window function 
may be applied at the spanwise boundaries as suggested by Guo et al. [8]. The window function renders both the 
mean flow and disturbances periodic along the z direction. Given a fundamental spanwise wave number, a single 
Fourier series may be used to represent both the mean flow and disturbances simultaneously. In this way, Eq. (10) 
reduces to a form similar to the 2D nonlinear PSE discussed in Chang [2], The major drawback of this approach is 
that a Fourier decomposition for the mean flow must be done with a proper choice of the fundamental spanwise 
wave number that is commensurate with the disturbance wave number of interest. The disturbance decomposition 
can only be constructed based on the same Fourier series. Even a linear problem would require solutions of the 
nonlinear PSE for this window PSE approach. Alternatively, the disturbance plane may be formulated in physical 
space; however, doing so would require more grid points on the (>’, z) plane to resolve the disturbance structure. The 
quasi-parallel eigenvalue solution at a given streamwise location may be used as the initial condition for the PSE 
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marching. However, it should be noted that such parallel initial condition would result in transient effects on the 
marching solutions. 

The plane-marching PSE approach is rather time-consuming despite the fact that it is a straightforward extension of 
the 2D PSE approach. Generation of the initial disturbance mode alone via a 2D eigenvalue solution alone would 
take significant computational time. It is therefore not suitable for routine engineering calculations for 3D boundary 
layers. Hu and Zhong [19] used the spectral collocation method to solve a form of governing equations similar to 
Eq. (10) for a 2D stability problem with periodic spanwise boundary conditions. Extensive studies on the plane- 
marching PSE would require more general boundary treatments on the cross plane. The plane-marching PSE method 
has not been implemented in LASTRAC at this moment. 

B. Surface-Marching PSE 

Instead of solving a plane full of disturbances simultaneously, in this approach, we look for a local PSE solution for 
a single surface-normal line. This local solution is made possible by evaluating x- and r-derivatives in Eq. (10) by 
using only neighboring surface points. By accounting for both derivatives on the body surface, we hope to include 
the effects of the mean flow and disturbance variation in X and z simultaneously. The concept is quite similar to the 
formulation for numerical solutions of the 3D boundary layer equations in which viscous terms in both x and z 
directions are dropped by orders of magnitude considerations and the governing equations are rendered parabolic 
along the major flow direction. The marching solutions can be performed along any direction on the body surface as 
long as the point under consideration falls within the zone of dependence formed by the known solutions given at 
neighboring points (see Wang [18]). The concept of zone of influence and dependence comes from characteristic 
theory for the simplified boundary layer equations. For a given point on the 3D surface, the zone of dependence is 
the wedge-shaped region formed by the intersecting inviscid and surface limiting streamlines. Typically, numerical 
solutions of the 3D boundary layer equations are obtained by marching on the 3D surface along both streamwise (x) 
and spanwise (r) directions (e.g., BL3D code [20]). 

To account for a more general disturbance on a 3D surface, we modify the disturbance form, Eq. (9), to 

~ iff a(£,z)d£+ f P{x 1 r/)dri\ 

V m =y' m (x.y t z)eV« K > (11) 


where a complex wave number f5 — f){x,z) is introduced in addition to a = a(x,z). The exponential part of 

Eq. (11) can be generalized to e w . The wave numbers then are defined by a ~ d0 / dx and P'—d6tdz, 
respectively. Mack [21] pointed out that within the framework of conservative wave theory, these two wave 
numbers are not entirely independent. To keep the wave front potential 9 intact, the following irrotational relation 
must be satisfied: 


da _ dp 
dz dx 


( 12 ) 


For 2D or infinite swept wing boundary layers, a constant and real spanwise wave number satisfies the above 
relation because the streamwise wave number is independent of z. 


The governing equation takes a form similar to Eq. (10), a set of nearly parabolic equations. Similar to the 2D PSE, 
the weak ellipticity associated with the pressure shape function gradient remains. A procedure similar to 3D 
boundary layer equations may be used to solve Eq. (10) over a rectangular domain on the surface provided proper 
initial conditions are given. Figure 2 depicts examples of the domain on a tapered wing surface. Unlike the case of a 
line-marching PSE method to be discussed later, the spanwise coordinate is fixed by the selection of the spanwise 
boundary line parallel to the leading edge. As shown in the figure, we need initial conditions on two surface lines: 
one parallel to the leading edge ( X = xo ) and one roughly aligned with the firee-stream direction ( z = zo ). The use 
of a nonorthogonal coordinate system offers a more flexible choice of initial condition lines. Both streamwise and 
spanwise wave numbers, a and p, are complex in general. Herbert [9] suggested that two normalization conditions 
derived from the kinetic energy norm can be used to determine these two wave numbers, a and (3. This approach has 
two drawbacks. Firstly, the irrotational relation Eq. (12) may not be satisfied. In addition, our numerical experiments 
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indicate that while this works for a weakly 3D boundary layer, two normalization conditions may compete with each 
other and result in divergence during the wave number iteration procedure for moderate to strong spanwise mean 
flow variations. 

Alternatively, we may use only the normalization condition in x to determine a and then impose Eq. (12) to evaluate 
p. This procedure implies that any change in the streamwise growth rate (determined by the imaginary part of a) 
along the spanwise direction will incur a non-zero imaginary part of p. For a moderately 3D boundary layer, the 
disturbance wave number (and its growth or decay) will change substantially during the PSE marching. This change 
results in numerical difficulties and leads to extremely large, unphysical growth or decay of the disturbance. 
Alternatively, we may require that only the real part of p be evaluated by Eq. (12), i.e., 

da, _ cfir 
dz ox 

The imaginary part of p is set to zero. This approach should not reduce numerical accuracy because the shape 
function dependency on z can account for any disturbance growth (or decay) and the irrotational wave front is 
preserved by the real part of p. Although not satisfying the wave front potential relation, a real constant P is also 
used in the 3D PSE marching procedure. As long as a large portion of the wave behavior is incorporated in p, 
additional phase or amplitude variation can be resolved by the shape function variation in z. Strictly speaking, the 
irrotational property of the wave vector, Eq. (12), is only valid for the wave form determined by the wave number 
norm chosen within the context of a nonparallel boundary layer. Additional shape function variations in both x and z 
would alter the wave front for different flow disturbances and for different wall-normal locations. Satisfying Eq. 
(12) is thus not a critical requirement for nonparallel boundary layers with strong mean flow variations in both 
surface directions. Further validations against DNS results are necessary to determine the merit of the different 
approaches discussed above. 

The surface marching procedure described above can account for disturbance evolution on the entire body surface in 
a consistent manner. However, it remains to be determined to what extent the marching solution is influenced by the 
initial conditions. In this research, we impose a local eigenvalue solution at (x 0 ,z 0 ) . A marching PSE solution is 

then obtained by neglecting the spanwise changes for the initial line that roughly aligns with the free stream (see 
Fig. 2). For the initial line nearly parallel to the leading edge, we impose the same eigenfunction obtained 
at (x 0 , z 0 ) . In general, this combination of initial conditions would result in transient effects for the marching 

solutions. However, the wave number iterations for a may damp out transient effects rapidly and limit the transient 
effects on the growth rate computed in the vicinity of the initial lines. The dependency of initial conditions that 
spans over a substantial part of the domain is the major drawback of the surface-marching PSE approach. Obtaining 
initial conditions from a higher-fidelity tool or advanced receptivity theories would help relieve this stringent 
requirement. Alternatively, a local line-marching PSE procedure with additional approximations can be used, which 
will be discussed in the following section. 

C. Local Line-Marching PSE 

In this approach, instead of solving for the disturbance evolution over a rectangular region on the surface, PSE 
marching is performed locally along a pre-identified line path on the surface. Along this line, a local PSE solution is 
solved for with the effects of spanwise mean flow variation being accounted for or neglected entirely. The main idea 
is to identify a path (e.g. streamline or group velocity line) that has relatively weak spanwise effects such that a local 
PSE solution represents most of the disturbance evolution. Governing equations and wave decomposition are similar 
to Eqs. (10) and (11), respectively. The spanwise wave number p is assumed to be real. By neglecting viscous 
second derivatives along z, coupled discretized PSE solutions at two neighboring spanwise wall-normal lines are 
solved simultaneously. This numerical technique is similar to the one used for the local nonparallel eigenvalue 
solution in a nonparallel boundary layer (see [2]). A real p does not pose any restriction because as in the surface- 
marching PSE, additional phase and amplitude modulations associated with the spanwise mean flow variation are 
captured in the shape function variation in z. Two separate norms in x are used to determine two streamwise wave 
numbers. Similar to the surface-marching case, a norm in z to determine a complex p results in a competing situation 
and wave number iterations would fail to converge. Mean flow variations along both directions on the surface are 
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accounted for in this procedure. The inviscid streamline is usually used as the integration path because the group 
velocity line requires a separate quasi-parallel linear stability solution. The spanwise (:) coordinate is aligned with a 
direction in which the mean flow variation is weak. The isobars seem to be a logical choice but it takes extra 
calculations to identify them. Numerical experiments indicate that the conical direction formed by the leading edge 
and trailing edge of a swept wing works fine for most cases. 

When a “good" spanwise direction is used, the local line-marching method may be further simplified by neglecting 
the - variation locally. This simplification is justified only as an engineering practice. Nonetheless, for many 
practical applications, this simplified local line-marching (without accounting for spanwise mean flow variations) 
offers a satisfactory solution that agrees well with those from the original local line-marching or the surface- 
marching approach. For example, in a tapered wing boundary layer, the conical ray direction may be used as the 
spanwise direction along which the mean flow variation is almost negligible. This s direction varies locally at each 
marching location and may be treated easily within the framework of a nonorthogonal coordinate system. For 
practical applications, the local line-marching method may be used for transition correlations due to its numerical 
efficiency. 

D. Quasi-parallel LST 

By further imposing the locally parallel assumption for the governing equations, an eigenvalue problem may be 
formed locally for a given disturbance frequency and complex spanwise wave number. In the saddle point method 
[3, 22], N-factor integration is performed along the group velocity direction obtained by solving the dispersion 
relation in conjunction with the requirement that da / dfi is real. To close the problem, it is further assumed that 

the disturbance growth rate is maximal along the group velocity direction. The spanwise wave number p is complex 
and varying along the marching direction in this case. However, the N-factor integration does not take into account 
the spanwise growth (or decay). Besides the saddle point method, the N-factor can be calculated by integrating the 
most unstable disturbance along the inviscid streamline. In many cases, there is little distinction in the N-factor 
contours between the saddle point and streamline maximization procedures. Both approaches are implemented in the 
LASTRAC.3d code. 


III. Mean Flow Computation 

Depending on the complexity of the 3D boundary layer configurations, various methods may be used to compute the 
mean flow. For a tapered wing, the conical assumption can be used to simplify the boundary layer equations (see 
[16]). Under this assumption, a similarity solution as a function of the conical ray angle and the wall-normal 
distance exists. Solutions at different spanwise locations can be constructed according to the conical radius 
measured from the origin of the rays formed by the leading edge and trailing edge to the location of interest. The 
boundary layer code BLSTA[16] is used to generate the transonic tapered wing case discussed in the next section. 

Boundary layer solutions for full 3D cases require a marching solution similar to the one discussed previously for 
the 3D surface-marching PSE procedure. Given the inviscid solution on the body surface, the boundary layer 
equations are solved consecutively via marching along both streamwise and spanwise directions on the body 
surface. The boundary layer code BL3D [20] is used to obtain 3D boundary layer solutions for the present study. 
The wall-normal velocity is neglected in the stability calculations because this velocity component as computed by 
BL3D has not been validated. Navier-Stokes solutions, of course, may be used as an alternative. 

IV. Results and Discussion 

The nonorthogonal coordinate system used in LASTRAC offers flexible choices of streamwise and spanwise 
directions. The N-factor integration in any local line-marching methods as discussed previously may be performed 
via marching along the chordwise (streamwise) surface grid line, the inviscid streamline, or the group velocity 
directions. In the grid line or streamline marching, a real the wave number p is used and disturbance growth along 
the spanwise direction is either accounted for in the shape function variation or completely neglected as an 
engineering approach. Three spanwise directions are available in LASTRAC: along the spanwise grid line, the 
conical direction, or perpendicular to the inviscid streamline. For a tapered wing, the first two options are naturally 
suitable because the spanwise mean flow variation along that direction is usually small. The last option is acceptable 
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downstream of the leading edge but would result in large errors near the leading edge where the streamline normal is 
almost perpendicular to the leading edge. For forebody configurations such as a circular or elliptic cone, the 
azimuthal direction is perhaps a good choice for p. An arbitrary combination of the marching and span wise 
directions may be employed for a line-marching PSE approach 

The line-marching PSE method provides a localized solution at the given marching path. The N-factor distribution is 
determined along the given path and any history effect due to spanwise variations is neglected. The surface- 
marching PSE method, on the other hand, can be used to obtain entire surface distribution of N-factors. Both line- 
marching and surface-marching PSE methods are implemented in the LASTRAC code. In the following discussion, 
we first validate the nonorthogonal coordinate system by comparing its results with existing validated orthogonal 
solutions. Then N-factor calculations using LST or various PSE methods are presented for several 3D boundary 
layer configurations. All calculations presented herein are performed by accounting for the body curvature along 
both chordwise and spanwise directions. 


A. Validation of the Nonorthogonal Coordinate System 

To validate the nonorthogonal coordinate system and governing equations, we calculate the stationary crossflow 
instability wave over an incompressible infinite swept wing boundary layer by using several combinations of 
marching and spanwise directions. As shown in Fig. 3, three marching directions are identified: along the chord 

normal ( x 0 ), along the inviscid streamline ( X, ), and along the free stream ( x 2 ). In the following calculations, the 
spanwise wavelength is fixed for two spanwise directions: parallel to the leading edge ( z 0 ) and along the direction 
normal to the inviscid streamline ( z, ). 

The quasi-3D solution computed by the earlier release of LASTRAC via marching along the chord normal direction 
with an orthogonal coordinate system (i.e. the (x 0 ,- 0 ) coordinate) is used as the baseline. It is considered the 

“correct" solution because the governing equations under an orthogonal coordinate are well-suited for an infinite 
swept wing. Numerical solutions using this orthogonal coordinate system have been validated against direct 
numerical simulations and excellent agreement was achieved [15]. Three additional combinations were 
investigated: (x 2 ,r 0 ), (Xj ,Z 0 ) , and (x, , z x ) . They are denoted by I, II, and III. respectively. Figure 4 shows the 

resulting nonparallel (PSE) N-factors based on disturbance kinetic energy for a spanwise wavelength of 10 mm 
(which is the most unstable stationary crossflow disturbance). We compare N-factors instead of growth rates 
because the former is an invariant under coordinate transformation while the latter is not. Both cases I and II give 
excellent agreement with the baseline solution. On the other hand, marching along the streamline with a normal 
spanwise coordinate (case III) gives an incorrect solution. This result is not surprising because case III differs from 
the others where the spanwise direction is consistent with the baseline solution. Although not shown here, quasi- 
parallel LST calculations using various coordinates produce identical results. 

The above comparisons indicate that a streamline or free-stream direction marching with a spanwise coordinate 
parallel to the leading edge gives a correct solution and the streamline orthogonal marching (case III) should not be 
used in a 3D boundary layer. We can also calculate the nonparallel N-factors by using the surface-marching PSE 
method. Figure 5 shows the resulting N-factor contours on the wing surface for a rectangular nonorthogonal domain. 
The spanwise direction is aligned with the spanwise domain boundary (parallel to the leading edge). The N-factor 
contours are parallel to the leading edge for the infinite swept-wing surface as expected. The N values along a 
chord-normal cut are also plotted in Fig. 4 and excellent agreement with the baseline solution is evident. 

B. Tapered Wing boundary Layer 

A tapered wing boundary layer provides a more realistic example for the study of 3D effects. Unlike an infinite 
swept wing, the mean flow in a tapered wing boundary layer varies along the spanwise direction. The mean flow is 
calculated by using the BLSTA code (Wie [16]) using the conical assumption, i.e., the mean flow does not vary 
along a ray originating from the point formed by the intersection of the leading and trailing edges. The mean flow is 
computed with an orthogonal coordinate system where marching is performed along a circular arc orthogonal to 
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both leading and trailing edges. Stability calculations could also be formulated by a similar conical assumption. 
However, we did not pursue this approach and instead use a more general nonorthogonal 3D framework in this 
study. 

A Mach 0.8 flow over a tapered wing was calculated. The leading and trailing edge sweep angles are 28° and 14°, 
respectively. The surface grid for mean flow calculations is shown in Fig. 6. Also shown in the figure are three 
different marching paths formed by tracing a gridline along the streamwise direction, along the inviscid streamline 
and along the group velocity direction. We perform local line-marching PSE calculations for the former two 
marching paths with a conical spanwise coordinate varying locally at each location of the tapered wing. Figure 7 
shows the corresponding linear PSE solutions obtained for the most amplified traveling crossflow mode with a 
frequency of 1.5 kHz and a spanwise wavelength of 12 mm. Evidently, both marching paths give identical solutions 
and may be used for tapered wing stability calculations. 

Quasi-parallel LST calculations were also performed by maximizing the disturbance growth rate and integrating N- 
factors along the inviscid streamline. Figure 8 shows the resulting N-factor along with those calculated by using the 
conventional saddle point method along the group velocity line. Both methods give almost identical N values except 
at small x/c. For a swept wing boundary layer, maximizing growth rates along the inviscid streamline tend to give 
very close N-factor results with the traditional saddle point method. 

Surface-marching PSE calculations were also performed for this configuration. The resulting N-factor distribution, 
as shown in Fig. 9(a), roughly scales according to the Reynolds number (i.e., the distance from the leading edge). 
Figure 9(b) depicts N-factor distributions at various spanwise locations. For comparison, the N-factors obtained by 
the local line-marching PSE solutions are also shown in the figure. The agreement between these two PSE marching 
solutions is very good. The small discrepancy is mainly due to the small difference in the initial conditions. 
Although not shown here, spanwise disturbance correction in the local line-marching PSE calculation was also 
performed for the above traveling crossflow mode. It was found that for the present case, the correction is negligible. 

The above results indicate that the local line-marching or surface-marching PSE method gives almost identical 
solutions for the test case. Any of these PSE methods may be used as an alternative to LST for instability wave and 
transition calculations for a tapered wing boundary layer in order to more accurately account for the nonparallel and 
surface curvature effects. 

C. Supersonic Finite Swept Wing 

The next test case is a Mach 2 flow over a 56-degree finite swept wing. The computed surface Cp distribution and 
wing geometry are shown in Fig. 10(a). Except near the leading edge, the pressure distribution indicates substantial 
3D variations. The computed crossflow Reynolds number contours are shown in Fig. 10(b). The maximum R cf is 

around 1500. The inboard side tends to have larger crossflow Reynolds numbers. The engulfed negative crossflow 
Reynolds number region in the middle is associated with the pressure variation in that region. The BL3D code reads 
in the inviscid solution and finds the attachment line location near the leading edge. An approximated similarity 
solution is then calculated at the attachment line. Using this solution as an initial condition, the boundary layer 
profiles are computed over the entire surface. 

We do not intend to perform exhaustive transition prediction calculations for this configuration. Instead, we will 
investigate 3D disturbance evolution for the most amplified traveling crossflow mode on the wing surface. From 
quasi-parallel LST calculations, the most amplified mode has a frequency of 6 kHz and a spanwise wavelength of 16 
mm. We first performed surface-marching PSE calculations. A local nonparallel eigenvalue calculation around one 
of the leading edge comers was performed for the most unstable mode by neglecting the spanwise mean flow 
variations. The computed streamwise wave number, a. and eigenfunctions were then imposed as initial conditions 
along a line parallel to the leading edge. The initial conditions along the inboard chordwise line were obtained by a 
local line-marching PSE calculation. While the chordwise initial conditions are a reasonable approximation, the ones 
along the spanwise direction are certainly not in light of the strong mean flow variation from inboard towards the 
outboard direction. Better initial conditions might be obtained by using local nonparallel eigensolutions at each 
location along the spanwise initial line. 
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We first assess the effects of initial conditions. Figures 1 1(a) and 1 1(b) show the N-factor contours evaluated by the 
peak chordwise velocity perturbations for two calculations using the local nonparallel eigensolutions obtained at 
inboard and outboard comers, respectively. Although both plots show similar trends, they differ in details. The 
difference in these results illustrates the dependency on initial conditions for the surface-marching PSE approach. 

In the above calculations, we imposed a constant spanwise wave number p. As discussed in the previous section, a 
complex spanwise wave number can also be imposed and updated during the surface-marching procedure according 
to the irrotational phase condition, Eq. (12). To maintain numerical stability, a real irrotational condition, Eq. (13) 
may also be used. The use of the additional constraint on the spanwise wave number implies that the modal wave 
number may vary downstream of the marching. This can be seen in Fig. 12(a)- 12(c) where the real parts of the 
effective spanwise wave number (after accounting for the shape function variation along the spanwise direction) 
contours evaluated at the peak chordwise velocity locations are plotted for real constant p, complex p using Eq. (12), 
and real p using Eq. (13), respectively. The additional constraint causes the effective wave number to vary 
significantly from the mid-chord region on and extend all the way to the trailing edge. It appears to introduce 
additional length scales to the disturbance field. The corresponding N-factor contours for the latter two cases shown 
in Fig. 12 are plotted in Figs. 13(a) and 13(b). Compared to Fig. 11(b) that uses the same initial conditions, these 
results reveal more complex structures and very large N values are observed within the “pockets” near the trailing 
edge. Although the appearance of multiple length scales in the disturbance evolution might look unphysical, further 
validations against DNS are necessary to judge the merits of imposing the irrotational condition. 

The local line-marching PSE method was also used to compute the above case. The same inboard comer initial 
conditions given for Fig. 11(a) were used for the local line-marching calculations. The spanwise mean flow and 
disturbance variations can be either accounted for by solving two wall-normal lines simultaneously or neglected as 
an approximation. The marching calculation was performed along the streamwise (chordwise) gridline direction 
with a spanwise coordinate aligned with conical rays. Figures 14(a)-14(c) depicts the N-factor contours evaluated 
by the total disturbance kinetic energy for the surface-marching, local line-marching with spanwise corrections, and 
local line-marching without spanwise corrections, respectively. Despite differences in details, the overall agreement 
of the three methods is acceptable. A close-up comparison reveals an N-factor difference of about 1 or 2. It must be 
cautioned that these comparisons are for a single disturbance mode. The local line-marching PSE method tends to 
give consistent trends with the surface-marching method. Because a transition correlation based on N-factor is 
empirical in nature, it remains to be answered that whether the surface-marching method offers a better prediction 
tool. 


D. Axisymmetric Cone at an Angle of Attack 

Besides swept wing boundary layers, axisymmetric or elliptic cones at an angle of attack represent another group of 
3D boundary layers. We investigate a Mach 3.5 flow over a 5-degree half-angle cone with a 4-degree angle of 
attack. Boundary layer transition on this configuration has previously been investigated experimentally by King 
[17]. The mean flow was computed by using the CFL3D code with a 256x33x177 grid along the streamwise, 
azimuthal, and wall-normal directions, respectively. Figure 15 depicts the surface mesh and computed inviscid 
streamlines. The maximum crossflow Reynolds number is about 1600. 

We first performed linear PSE calculations along two different paths as shown in Fig. 16. The gridline path is 
roughly aligned with the ray of the cone generator. A spanwise wavelength of 1 mm (which corresponds to the most 
amplified stationary crossflow mode) was fixed along the azimuthal direction for both cases. Figure 17 shows the 
resulting nonparallel stationary crossflow instability N-factors (computed by a local line-marching PSE without 
spanwise corrections) measured by the disturbance kinetic energy. The large N values indicate that the flow is very 
unstable. Unlike the previous swept wing examples, two different marching paths give a difference in N-factor of 
about 5. A more comprehensive study is necessary to judge the superiority of either marching path. Previous results 
seem to indicate that marching along the inviscid streamline in general gives better transition correlation while 
marching along the ray usually overpredicts the N value. 

The above PSE results were obtained by neglecting the spanwise (azimuthal) mean flow and mode shape variation. 
Figure 18 shows linear PSE results obtained by streamline marching with and without accounting for the spanwise 
variation. As shown in the figure, the 3D effects result in a smaller N-factor and the difference is slightly larger than 


12 

American Institute of Aeronautics and Astronautics 



1. In the linear PSE calculation with spanwise correction, the effective spanwise wave number becomes complex 
because the shape function is allowed to vary in the spanwise direction. The maximum spanwise wavelength 
correction is only about 2 percent. It appears that the disturbance is decaying along the spanwise direction (from the 
windward pointing towards the leeward direction), which is expected because of the decrease of crossflow strength 
along that direction. Due to the lack of resolution along the azimuthal direction and large mean flow gradient near 
the apex of the cone, the surface-marching PSE method gives solutions that possess strong and unphysical transient 
effects near the leading edge. A better resolved mean flow and more carefully treated initial conditions are necessary 
for further studies. 

We also performed transition correlation for stationary crossflow instability by maximizing the growth rate with 
respect to spanwise wave number by using the conventional LST approach. The streamline marching path was used 
for integration. Figure 19 shows the location where the N value first reaches 10 for each streamline. The results 
show that transition location moves upstream near the center plane. Both windward and leeward planes are very 
stable to stationary crossflow disturbance. These results seem to be consistent with previous investigations (Malik 
and Balakumar [5]). 


V. Conclusion 

Advanced transition prediction methods based on traditional and the state-of-the-art PSE methods are studied for 3D 
boundary layers. Numerical formulations of various transition prediction methods for general 3D boundary layers 
are discussed. The surface-marching and local line-marching PSE approaches provide a means of calculating 
disturbance evolution in 3D boundaiy layers beyond the traditional quasi-parallel LST methods. The use of a 
nonorthogonal coordinate system facilitates easy adaptation to the mean flow grid and flexibility in choosing the 
combination of spanwise and marching directions. These 3D transition prediction methods are implemented in the 
new release of the LASTRAC code. 

Several test cases that cover a range of speed regimes are studied by using the surface-marching and line-marching 
PSE methods. The results indicate that both methods can capture 3D effects due to mean flow variations. For a 
tapered wing boundary layer, the line-marching PSE with a conical spanwise direction gives identical solutions to 
that by the surface-marching method. For the supersonic swept wing case, initial conditions have noticeable 
influence on the computed N-factor contours. However, with the same initial conditions, line-marching and surface- 
marching PSE methods give qualitatively similar but quantitatively different solutions in N-factors. Further 
validations against existing data or DNS results are necessary to assess the merits of each approach. 
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Figure 1: Nonorthogonal body-fitted 

curvilinear coordinate over a 3D surface. 


Figure 2: Rectangular domains for 

surface-marching PSE calculations. 





Figure 3: Three different marching 

coordinates for an infinite swept wing 
boundary layer. 



Figure 4: Comparison of PSE N-factors for a 
stationary crossflow mode using various 
coordinate systems and surface-marching PSE 
for an infinite swept wing boundary layer. 
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Figure 5: N-factor distribution calculated by 
surface-marching PSE for the infinite swept- 
wing boundary layer using a nonorthogonal 
domain. 



Figure 6: Geometry and surface grid for a 
transonic tapered wing, showing three 
marching paths and the conical spanwise 
coordinate. 



Figure 7: Comparison of linear PSE N-factors 
for /= 1 .5 kHz and X = 12 mm using two 
marching paths for a transonic tapered wing 
boundary layer. 



Figure 8: Comparison of maximized LST N- 
factors for /= 1 .5 kHz using two integration paths 
for a transonic tapered wing boundary layer. 
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Figure 9: Nonparallei N-factor (based on disturbance kinetic energy) distribution for 
/= 1-5 kHz. X = 12 mm on the transonic tapered wing surface obtained by surface- 
marching PSE: (a) N-factor contours on the surface; (b) N-factors from local line- 
marching at various spanwise locations. 


Y 



(a) 


(b) 


Figure 10: Pressure coefficient (a) and crossflow Reynolds number (b) distributions on the 
surface of a Mach 2, 56-degree finite swept wing 
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Figure 1 1: Surface-marching N-factor contours for/= 6 kHz and X = 16 mm using two 
different initial conditions: (a) inboard corner: (b) outboard corner. 



Figure 12: Real part of the effective spanwise wave number evaluated at peak chordwise 
velocity locations using: (a) constant real P; (b) Eq. (12); (c) Eq. (13). 
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Figure 13: Surface-marching N-factor evaluated at peak chordwise velocity locations 
using: (a) Eq. (12); (b) Eq. ( i 3 ). 
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Figure 14: Comparison of N-factor evaluated by total kinetic energy obtained by (a) 
surface-marching: (b) local line-marching with spanwise correction: (c) local line- 
marching without spanw ise correction. 
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Figure 15: Surface mesh and streamlines for a Figure 16: Two marching paths on a Mach 3.5 

Mach 3.5 flow over a cone at 4-degree angle of flow over a cone at 4-degree angle of attack, 

attack. 



Figure 17: Linear PSE N-factors 
computed by two different marching 
paths on the Mach 3.5 cone surface. 



Figure 18: N-factor calculated by linear PSE 
using local line-marching along the inviscid 
streamline with and without spanwise mean 
flow and disturbance correction. 
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N = 10, Stationary Crossflow 


Figure 19: Stationary crossflow N-factor distribution 
calculated by LST marching along the streamline with 
maximizing the local disturbance growth rates. 
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